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Nomenclature 


ARle 

Aspect ratio at leading edge 

k 

Mesh gradation coefficient 

Cf,x 

Skin friction coefficient component 

r 

Distance from reentrant comer tip 

C p 

Pressure coefficient 

P 

Convergence order 

Cl 

Lift coefficient 

Q 

Exact solution of the Laplace equation 

Cm 

Pitching moment coefficient 

s 

Size of reentrant-corner cell 

Cd 

Total drag coefficient 

u , w 

Velocity components 

Cd p 

Pressure drag coefficient 

X, z 

Coordinates relative to reentrant corner tip 

C Dv 

Viscous drag coefficient 

x,y,z 

Orthogonal coordinates in streamwise, 

M 0 o 

Freestream Mach number 


spanwise, and normal directions 

N 

Number of degrees of freedom 

XTE.ZTE 

Reentrant corner / trailing-edge position 

Pref 

Reference pressure 

z + 

Nondimensional boundary layer spacing 

Pt 

Total pressure 

a 

Angle of attack 

Pr 

Prandtl number 

a 

Inverse relative exterior angle 

Pr t 

Turbulent Prandtl number 

p 

Stretching factor 

Re 

Reynolds number 

e 

Discretization error 

Pref 

Reference temperature 

0 

Polar coordinate angle 

T t 

Total temperature 

K 

MUSCL scheme parameter 

a re f 

Reference speed of sound 

Pref 

Reference eddy viscosity 

c 

Chord length 

Pt 

Eddy viscosity 

havg 

Averaged mesh spacing 

Voo 

Freestream kinematic viscosity 

h > heff 

Characteristic mesh size 


Mapping coordinate 

hmax 

Maximum mesh size 

UJ 

External angle for reentrant corner 


I. Introduction 

With ever-increasing computing power and recent advancements in solver technology, turbulent flows are routinely 
simulated on high-density grids with many millions of degrees of freedom. While accurate and reliable Computational 
Fluid Dynamics (CFD) analysis of attached turbulent flows is now possible, accuracy and robustness of separated 
turbulent flow simulations in complex geometries are still not adequate. Errors in such simulations are typically 
attributed to three sources. (1) The modeling error is due to approximations in the continuous formulation (e.g., 
in differential equations describing turbulent flows, or in geometry definitions, or in boundary conditions) and is 
defined as the difference between the exact continuous solution of the model formulation and the real flow. (2) The 
discretization error is due to approximations in discretizing the continuous formulation on a specific grid and is defined 
as the difference between the exact discrete and continuous solutions. (3) The iterative ( algebraic ) error is due to an 
imperfect iterative solution process for the discrete formulation and is defined as the difference between the exact and 
approximate solutions to the discrete formulation. 

A Turbulence Modeling Resource (TMR) website 1,2 has been recently established at NASA Langley Research 
Center to describe, standardize, verify, and validate formulations of common turbulence models. The purpose of this 
website is to avoid ambiguity associated with specific implementations of turbulence models in CFD codes. The 
turbulent flows considered in this paper are modeled by Reynolds-averaged Navier-Stokes (RANS) equations and the 
one-equation linear eddy-viscosity Spalart-Allmaras (SA) turbulence model. 2 3 The formulation details are available at 
the TMR website. 

This paper reports on an attempt to eliminate (or at least minimize) the discretization and iterative errors by con- 
ducting an extensive grid convergence study for relatively simple benchmark turbulent flows. Current guidelines for 
grid convergence studies 4,5 emphasize a parametric similarity of grids forming a family and an asymptotic convergence 
order, which is expected to be observed on three-to-four fine grids in a family. The mesh resolution required for estab- 
lishing a convergence order is sought through a uniform grid refinement. For structured grids, a family of uniformly 
refined grids is typically derived recursively, starting from the finest grid. Each coarser grid in the family is derived 
from the preceding finer grid by removing every-other grid plane/line in each dimension. This uniform-refinement 
approach has rigorous mathematical foundations. However, it is also expensive because it lacks the flexibility of a 
local refinement, which is the basis for effective grid adaptation methods. 

Ideally, solutions obtained on grids in a family would monotonically converge in grid refinement to the continuous 
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solution with the same asymptotic convergence order for any solution quantity. Then, the entire solution can be 
extrapolated to the limit of the infinite grid refinement. 

A single asymptotic convergence order characterizing the entire solution is typically observed only for model 
problems with smooth continuous solutions. In practical computations, an asymptotic convergence order remains 
elusive. Computations presented at the recent Drag Prediction Workshop (DPW-V) 6,7 illustrate the difficulty in real- 
izing asymptotic convergence order for turbulent flows. A relatively close agreement was observed between different 
computations in predicting the drag on a three-dimensional (3D) wing-body configuration, but the grid convergence 
properties were very different between solutions computed with different codes or even between same-code solutions 
computed on different grid families. In anticipation of second-order convergence, the values of drag were plotted 
versus N~ 3 , where N is the number of degrees of freedom, but second-order convergence was not observed. 

The study reported in this paper aims at computing highly accurate reference solutions for some benchmark 
turbulent-flow cases and at providing some guidance on accuracy variation for grid families with different mesh den- 
sity patterns. The main grid-convergence test case considered is a turbulent flow around the 2D NACA 0012 airfoil 
at 10° angle of attack. The RANS equations are solved on uniformly -refined, structured, high-density grids by three 
well-established CFD codes: CFL3D (NASA), FUN3D (NASA), and TAU (DLR), which use different discretization 
and iteration schemes. Advanced turbulent-flow solver technologies recently implemented in two of these codes 8-10 
provide means for minimizing effects of iterative errors. FUN3D and TAU converge all residuals, including the resid- 
ual of the SA turbulence-model equation, to levels comparable with the machine zero. CFL3D converges the density 
residual to the level of 10“ 13 and the SA model residual to the level of 10“ 7 ; the corresponding aerodynamic forces 
converge to at least five significant digits; and the pitching moment converges to at least four significant digits. 

The study began with an attempt to characterize the grid-refined solutions by using grids offered in the “Turbu- 
lence Model Validation Cases and Grids” section of the TMR website. The grid-refinement study was conducted 
using FUN3D solutions on the family of grids then available 9 — the finest grid had about 1M degrees of freedom. 
Convergence of the lift and pitching moment coefficients did not exhibit any clear order property. The lift values in- 
creased with grid refinement and then decreased on the finest grid. The pitching moment coefficients were continually 
increasing with grid refinement, but showed no consistent order. A detailed inspection of the solutions on the surface 
of the airfoil revealed erratic convergence of the pressure coefficients near the trailing edge. 

This observation motivated the current grid convergence study that involves solutions obtained with the three CFD 
codes on three expanded grid families. The finest grids in each family have about 15M degrees of freedom. The grid 
families differ in the trailing-edge resolution and are now available in the “Cases and Grids for Turbulence Model 
Numerical Analysis” section of the TMR website. Convergence sensitivities to the trailing-edge resolutions as well 
as to various discretization aspects, such as grid elements and the order of approximation for the turbulence-model 
convection term, have been considered. 

Besides the NACA 0012 study, an existing flat plate test case has been extended and used to study grid convergence. 
FUN3D solutions have been computed on a set of structured grids. The grids are also available in the “Cases and 
Grids for Turbulence Model Numerical Analysis” section of the TMR website. Methodologically, the paper follows 
the current guidelines for grid convergence study: families of uniformly refined grids are used and convergence orders 
of local and global solution quantities are reported. 

Additionally, the convergence degradation long-known for solutions of elliptic (pure diffusion) equations on uni- 
formly refined grids with geometric singularities 11 is revisited. Elliptic equations describe diffusion phenomena and 
thus apply directly to low-Reynolds number (i.e., Stokes) flows. Near surfaces, turbulent- flow solutions are expected 
to be similar to Stokes-flow solutions. For a sharp trailing edge, the pure-diffusion solution exhibits a square-root 
behavior near the trailing-edge singularity and has unbounded derivatives. The discretization error converges on uni- 
formly refined grids with the first order in the L\ norm and with an order of 0.5 in the L ^ norm. A series of structured 
(non-uniformly) refined grids that have a higher refinement rate near the singularity than in the rest of the domain 
are developed and shown to recover the convergence rate obtained for smooth flows on domains without singularities. 
This grid refinement strategy has not been applied to turbulent-flow computations. However, its success in applica- 
tion to the pure-diffusion equation provides an indication that improved resolution near geometric singularities would 
improve the accuracy of turbulent-flow simulations. 

The material in the paper is presented in the following order. First, the CFL3D, FUN3D, and TAU codes used 
in the study are described in Section II, including discretization details and iterative convergence strategies. Then, a 
benchmark turbulent flow around the NACA 0012 airfoil is introduced in Section III together with a description of the 
grid families and the corresponding numerical solutions obtained with the three codes. A detailed description of the 
solutions and grid convergence is provided for future verification of CFD solvers. A study of solution sensitivity to 
the variation of discretization methods and grid elements is presented. Section IV reports on a grid convergence study 
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for a turbulent flow around a finite flat-plate configuration. Finally, concluding remarks are offered. The appendix 
considers effects of geometric singularity on accuracy for solutions of a pure-diffusion equation and presents some 
strategies to overcome the convergence degradation. 

II. CFD Codes Used in the Study 

This section introduces the three well-established large-scale practical CFD codes used in the grid-convergence 
study for the NACA 0012 airfoil. The codes represent the state of the art in aerodynamic computations. 

A. CFL3D 

CFL3D is a structured-grid multi-block cell-centered finite- volume code widely applied for analysis of complex flows. 
It has been used in many recent workshops involving complex turbulent flows 12-15 and for computing benchmark 
turbulent-flow solutions at the TMR website. It uses second-order, upwind-biased spatial differencing scheme (a 
MUSCL scheme 16,17 corresponding to k = 1/3 that allows a third-order accuracy in one dimension (ID)) for the 
convective and pressure terms, and second-order differencing for the viscous terms; it is globally second-order ac- 
curate. Roe’s flux difference- splitting method is used to obtain fluxes at the cell faces. The option to model the full 
Navier-Stokes meanflow equations is exercised for all cases. In distinction from the other two codes that use the 
SA-neg scheme, 18 CFL3D uses the standard SA one-equation turbulence model for this study. The negative values 
of the Spalart turbulence variable are not allowed; the minimum values are clipped at 10 -12 . First and second-order 
approximations for the convection term in the turbulence-model equation are available. A second-order approximation 
is used for the NACA 0012 case on grids of Family II (see Section III for grid family definitions). Initially, the first- 
order approximation was used for all computations; however, the second-order approximation was found to make a 
significant difference on Family II grids (see figures on the TMR website). In this study, the first-order approximation 
is used for computations on grids of Family I and Family III. The turbulence-model diffusion term uses the thin-layer 
approximation. The iteration scheme is loosely coupled, i.e., first, the meanflow equations are advanced with the eddy- 
viscosity fixed, then the turbulence-model equation is advanced with the meanflow solution fixed. CFL3D employs 
local time-step scaling, grid sequencing, and multigrid to accelerate convergence to steady state. 

B. FUN3D 

FUN3D is a finite- volume, node-centered, unstructured-grid RANS solver, which is also widely used for high-fidelity 
analysis and adjoint-based design of complex turbulent flows. 15,19-25 FUN3D solves governing flow equations on 
mixed-element grids; the elements are tetrahedra, pyramids, prisms, and hexahedra. At median-dual control- volume 
faces, the inviscid fluxes are computed using an approximate Riemann solver. Roe’s flux difference splitting is used 
in the current study. For second-order accuracy, face values are obtained by a MUSCL scheme, with unweighted 
least-squares gradients computed at the nodes. If grid lines are available, e.g., within boundary layers or in the wake, 
there is an option to use a directional gradient that exploits a ID line-mapping along the grid lines. For this study, the 
MUSCL scheme coefficient is set to k = 0.5 for the meanflow equations and to k = 0 for the turbulence convection 
term. The viscous fluxes use full viscous stresses. For tetrahedral meshes, the viscous fluxes are discretized using the 
Green-Gauss cell-based gradients; this is equivalent to a Galerkin type approximation. For non-tetrahedral meshes, the 
edge-based gradients are combined with Green-Gauss gradients; this improves the h-ellipticity of the viscous operator. 
The diffusion term in the turbulence model is handled in the same fashion as the meanflow viscous terms. FUN3D 
uses the SA-neg variant of the SA turbulence model 18 that admits negative values for the Spalart turbulence variable. 
This variant was designed for improved numerical behavior. The SA-neg model is identical to the original SA model 
for positive values of the Spalart turbulence variable. FUN3D uses a second-order approximation for the convection 
term in the turbulence-model equation. 

A multigrid solver is used to converge residuals. The relaxation scheme in this multigrid solver is a hierarchical 
nonlinear scheme. On the innermost level it uses a preconditioner based on a defect-correction method and iterates on 
a simplified first-order Jacobian with a pseudo-time term. One preconditioner iteration involves an implicit-line pass 
through the portion of the domain where implicit grid lines are defined, followed by a point-implicit sweep through 
the entire domain. The number of preconditioner iterations may vary for different nonlinear iterations. This variable 
preconditioner is wrapped with a Generalized Conjugate Residual (GCR) method to form a Jacobian-free linear solver 
that uses Frechet derivatives to approximate linear residuals. A nonlinear controller assesses the correction computed 
by the linear solver. The controller is responsible for the CFL adaptation strategy and for deciding when to update 
the Jacobian. As a result of this assessment, the suggested correction can be applied fully, partially, or completely 
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discarded; the current Jacobian may be updated or reused in the next iteration; and the current CFL number may 
increase, decrease, or stay the same. In the relaxation scheme, the iterations can be tightly or loosely coupled, i.e., 
applied to the meanflow and turbulence equations collectively or separately. The multigrid iterations are always 
coupled in the sense that the meanflow and turbulent equations are solved on coarse grids and the meanflow and 
turbulence variables use a coarse-grid correction. Initially, the CFL number is ramped over a prescribed number of 
iterations, but then it automatically changes within prescribed bounds. The coarse-grid corrections are also assessed 
by the nonlinear controller and can be applied fully, partially, or completely discarded. 

C. TAU 

TAU is a finite- volume node-centered unstructured-grid RANS solver widely used for a broad range of aerodynamic 
and aero-thermodynamic problems. 26 It offers coupling interfaces to other disciplines like structure and flight mechan- 
ics to allow for multidisciplinary simulations. 27 A full derivative is available for adjoint-based shape optimization. 
TAU solves the 3D compressible time-accurate RANS equations on grids with mixed elements, including tetrahedra, 
pyramids, prisms, and hexahedra. Control volumes are constructed by median-dual partition. The numerical scheme 
is based on a second-order, finite- volume formulation. At control volume faces, the inviscid fluxes are computed using 
a central difference scheme with an added matrix- valued artificial viscosity. 10 To deal with highly stretched meshes, 
a cell stretching coefficient is included into the scheme. The full viscous fluxes of the meanflow and turbulence equa- 
tions are discretized using an edge-normal gradient formulation as an augmented average of the adjacent Green-Gauss 
cell gradients. 28 Various turbulence models are available, ranging from eddy viscosity to full differential Reynolds 
stress models, 29 including options for Large-Eddy Simulation (LES) and hybrid RANS/LES. The SA-neg model 18 
is used as the turbulence model in this study, and the SA model convection term is discretized using a second-order 
approximation. 

A multigrid solver based on agglomerated coarse grids is used to converge to steady state. The baseline relaxation 
scheme of TAU in this multigrid solver is an implicit Lower-Upper Symmetric Gauss-Seidel (LU-SGS) scheme. 30 
Recently the LU-SGS scheme was embedded in a hierarchy of smoothers derived from a general implicit Runge- 
Kutta method to further improve reliability and efficiency of the solutions algorithms of TAU. 10 The smoothers, e.g. 
a first-order preconditioned Runge-Kutta or Newton-Krylov generalized minimal residual (GMRES) methods, can 
be considered as simplified Newton methods. The smoothers differ in Jacobian approximations and in the solution 
methods used for the arising linear systems. 

III. NACA 0012 Airfoil 

A grid convergence study for a turbulent flow around the NACA 0012 airfoil is presented in this section. This test 
case corresponds to the NACA 0012 case in the “Cases and Grids for Turbulence Model Numerical Analysis” section 
of the TMR website. The goals of this study are (1) to establish an accurate reference solution that can be used for 
verification of CFD solutions computed with the SA turbulence model, (2) to evaluate the effects of grid resolution 
near a sharp trailing edge on convergence of turbulent-flow solutions, and (3) to assess sensitivity of the solutions to 
variations of discretization methods and grid elements. 

A. Flow Parameters, Boundary Conditions, and Discretization Details 

A turbulent essentially incompressible (M 0 0 =0.15) flow around the NACA 0012 airfoil at a = 10° angle of attack is 
considered. For the purposes of this study, the definition of the NACA 0012 airfoil is slightly altered from the original 
definition, so that the airfoil closes at c = 1 with a sharp trailing edge. The exact definition is available at the TMR 
website. 1 The Reynolds number computed per chord length is Re = 6M. The computational domain and boundary 
conditions are sketched in Fig. 1. The farfield boundary conditions are based on inviscid characteristic methods. A no- 
slip adiabatic wall condition is specified at the airfoil surface. FUN3D and TAU have a strong implementation of the 
wall boundary conditions; CFL3D has a weak implementation of the wall boundary conditions. T re f = 540° R is the 
freestream static temperature. Although a farfield point vortex boundary condition correction 31 is recommended at the 
TMR website, the results below are presented without such a correction. This simplification facilitates comparisons 
with emerging high-order and mesh adaptation capabilities. 32,33 The farfield value of the Spalart turbulence variable 
is v far field = The Prandtl number is taken to be constant at Pr = 0.72, and the turbulent Prandtl number is 

taken to be constant at Pr t = 0.9. The molecular viscosity is computed using Sutherland’s Law. 34 
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NACA 0012 Boundary Conditions, 



Figure 1. Domain and boundary conditions. 


Table 1. NACA 0012: summary of along-surface mesh spacing on the finest 7, 169 x 2, 049 grids in three families. 


Grids x « 0 

(Leading-edge) 
Family I 0.0000125c 

Family II 0.0000125c 

Family III 0.0000125c 


x « 1 

(Trailing-edge) 

0.0001250c 

0.0000125c 

0.0000375c 


x « 0.5 

(Middle of the surface) 
0.00123c 
0.00155c 
0.00139c 


B. Grids 


449 x 129 grid, far view 



-500 0 500 


449 x 129 grid, near view 



(a) Far view. 


(b) Near view. 


Figure 2. Computational domain and a 449 x 129 grid for NACA 0012 airfoil. 

Three families of grids are generated with a farfield extent of approximately 500 chord lengths. Figures 2 (a) 
and (b) show two views of the 449 x 129 grid of Family I. Family I grids have the density distribution similar to the 
distribution used on grids of the family available on the TMR website prior to this study. Family II grids are clustered 
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(a) Family I. 


(b) Family II. 


(c) Family III. 


Figure 3. Near view of trailing-edge grids of Families I, II, and III. 


near the trailing edge and Family III grids are intermediate between the two. See Figs. 3 (a), (b), and (c) for near- 
trailing-edge views of the 449 x 129 grids for each family. A series of seven nested grids are generated for each family, 
ranging from the finest 7, 169 x 2, 049 to the coarsest 113 x 33 grid. The grid topology is a so-called “C-grid.” Each 
of the grids wraps around the airfoil from the downstream far field, around the lower surface to the upper, then back to 
the downstream far field; the grid connects to itself in a one-to-one fashion in the wake. There are 4097 points on the 
airfoil surface on the finest grid, with 1537 points along the wake from the airfoil trailing edge to the outflow boundary. 
Each family’s finest grid has the minimum normal spacing at the wall of 10“ 7 . The spacing along the airfoil surface 
is documented in Table 1 . The leading-edge spacing is the same for all families and corresponds to the aspect ratio of 
ARle = 125. The trailing-edge spacing is largest for the Family I grids and ten times larger than the leading-edge 
spacing. On Family II grids, the trailing-edge spacing is the same as the leading-edge spacing. On Family III grids, the 
trailing-edge spacing is between that of Family I and Family II and three times larger than the trailing-edge spacing of 
the corresponding Family II grids. The family name convention is not consistent with the variation of the trailing-edge 
mesh spacing. The authors choose to keep the same names for grid families as in the “Cases and Grids for Turbulence 
Model Numerical Analysis” section of the TMR website. 

The mesh spacing in the middle of the airfoil surface changes between the families. The along-surface spacings at 
x « 0.5 are 0.00123c, 0.00139c, and 0.00155c for families I, III, and II, respectively. The corresponding aspect ratios 
are 12300, 13900, and 15500. The relative increase in the middle-of-chord mesh spacing and aspect ratio between 
families I and II is approximately 25%. The middle-of-chord aspect ratios are approximately two orders of magnitude 
higher than those at the leading edge. 

C. Grid Convergence of Aerodynamic Coefficients 

This section reports on convergence of aerodynamic coefficients on grids of families I, II, and III. Figures 4-8 compare 
convergence of the total drag ( Cd ), pressure drag ( Cd p ), viscous drag ( Cd v ), lift (Cl), and pitching moment (Cm) 
with respect to the quarter-chord location. The variations are shown versus a characteristic mesh spacing h = y/l/N. 
FUN3D computations are shown only on the four finest grids of Family III. To accommodate a detailed scale for the 
Cl and Cm coefficients, only results on the three finest grids in each family are shown in Figs. 7-8. 

Convergence plots of drag coefficients shown in Figs. 4-6 are similar on grids of different families. Convergence 
plots of lift and pitching moment differ significantly between grid families. Note, however, that, relatively speaking, 
the vertical scale for the lift figures is significantly smaller (showing variations in the fourth significant digit) than 
vertical scales for the drag and moment figures (showing variations in the third and first significant digits, respectively). 
Relatively large deviations of the CFL3D lift and moment coefficients from the corresponding FUN3D and TAU 
coefficients observed on Family I and Family III grids are partially explained by variations in the discretization scheme 
for the SA model equation. Recall that CFL3D solutions are computed with the first-order approximation for the 
convection term in the SA model equation on grids of Family I and Family III. Although not shown here, results on 
the TMR website demonstrate the effect of the SA model discretization order when using Family II grids. 

All aerodynamic coefficients are predicted with a small variation between all the three codes on the finest grids of 
all families: the drag variation between the codes and grid families is less than 1 count (less than 1%), the lift variation 
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is less than 0.2%, and the pitching moment variation is less than 10%. The maximum variation is observed between 
the three codes on the Family I grid; the corresponding variation on the Family II grid is an order of magnitude smaller: 
0.1 count (less than 0.1%) for Cd, 0.005% for Cl , and 1.5% for Cm- 

The variations between the codes on finer grids within the same family are also smaller on the Family II grids than 
on grids of other two families. In fact, the lift and pitching moment coefficients computed on the Family I grids appear 
to be converging to values different from those computed on the Family II grids. This discrepancy motivated the 
introduction of an intermediate Family III. Family III solutions on coarser grids appear to be converging to yet another 
limit, but on the finest grid turn toward the values computed on the Family II grids. This behavior is observed for all 
codes and indicates that the solution variations due to differences in the trailing-edge resolution are larger than the 
variations due to differences in discretization schemes. Green dotted lines in Figs. 7-8 show the values corresponding 
to the infinite grid refinement computed by a linear extrapolation fitting the two finest grids. On grids of Family I, 
the extrapolated lift coefficients vary between values of 1.0885 and 1.0905 and the extrapolated pitching moment 
coefficients vary between 0.0069 and 0.0074. On grids of Family II, the extrapolated lift coefficient is 1.0910 and 
the extrapolated pitching moment coefficient is 0.0068. Note that the lift and pitching moment coefficients computed 
from presumably the most accurate solutions on the finest Family II grid lie outside of the range spanned by the lift 
and pitching moment values extrapolated from solutions on grids of Family I and Family III. 





(a) Family I. (b) Family II. (c) Family III. 


Figure 4. NACA 0012: Grid convergence of the total drag coefficient (Cd)- 





(a) Family I. 


(b) Family II. 


(c) Family III. 


Figure 5. NACA 0012: Grid convergence of the pressure drag coefficient ( Cd p )• 

Figures 9-11 show variations of forces and moment computed on grids of Family II with respect to h 2 = 1/N. 
The results are shown for forces and moment computed over the full airfoil and over the areas near the trailing and 
leading edges. The local integration areas are defined in Table 2. The right end of the leading-edge integration interval 
is selected as the ^-coordinate of the surface node on the 897 x 257 grid nearest to x = 0.1. Analogously, the left 
end of the trailing-edge integration interval is selected as the x-coordinate of the surface node on the same 897 x 257 
grid nearest to x = 0.9. These end nodes are present on four finer grids. The contributions to the forces and moment 
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(a) Family I. 


(b) Family II. 


(c) Family III. 


Figure 6. NACA 0012: Grid convergence of the viscous drag coefficient ( Cd v )• 





(a) Family I. 


(b) Family II. 


(c) Family III. 


Figure 7. NACA 0012: Grid convergence of the lift coefficient (Cl). 





(a) Family I. 


(b) Family II. 


(c) Family III. 


Figure 8. NACA 0012: Grid convergence of the pitching moment coefficient (Cm)* 


are much larger in the leading-edge region than in the trailing-edge region. The results for CFL3D computations on 
quadrilateral grids and for FUN3D computations on triangular grids are plotted only for the full airfoil. 

The convergence plots for lift, moment, and pressure drag are almost linear over the three finest quadrilateral 
grids for all three codes, indicating apparent second-order convergence. Lift and moment computed by FUN3D on 
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Table 2. Leading and trailing-edge integration areas. 


Leading Edge Trailing Edge 

0.0 < x < 0.100177952877727 0.899166466843597 < x < 1.0 


triangular grids over the full airfoil show apparent second-order convergence over the three coarser grids and a sharp 
turn on the finest grid. The values of Cl and Cm computed by FUN3D on the finest triangular grid are close to the 
values computed by FUN3D and TAU on the finest quadrilateral grid. The Cd p coefficient computed by FUN3D on 
triangular grids converges with an apparent order higher than second. Convergence plots for the viscous drag show 
less than second-order convergence for the drag computed by FUN3D and TAU over the full airfoil and over the 
leading-edge area. Variations of drag in the trailing-edge area appear very small. The extrapolated, grid-refined values 
of aerodynamic coefficients computed with different codes are not the same. CFL3D extrapolates lift and pitching 
moment to values somewhat different from the values extrapolated by FUN3D and TAU. These discrepancies may be 
a result of differences in implementation of the SA turbulence model. CFL3D employs a thin-layer approximation 
for the diffusion term and a standard SA formulation that does not allow negative values for the turbulence variables; 
FUN3D and TAU use a full-diffusion approximation and the SA-neg variant of the SA model. The extrapolated values 
of the lift and moment in the trailing-edge area show some differences between FUN3D and TAU solutions as well. 

D. Surface Pressure and Skin Friction 

This section compares the surface pressure and skin friction coefficients from the FUN3D, CFL3D, and TAU solutions 
on the finest 7, 169 x 2, 049 grid of Family II. In moderately zoomed views focused on the leading and trailing edges 
(Fig. 12), the solutions are indistinguishable. Only with a super zoom (Fig. 13) do some differences come to light. 
Figures 13 (a) and (b) compare solutions close to the minimum pressure and the maximum skin friction locations near 
the leading edge. The CFL3D solution shows a smaller pressure and less skin friction than the FUN3D and TAU 
solutions, which are indistinguishably close to each other, even on the super-zoom view. The largest differences are 
observed in the immediate vicinity of the trailing edge (Figs. 13 (c) and (d)); FUN3D and TAU solutions indicate a 
small area of a positive load, while the CFL3D solution indicates a negative load in this area. The TAU solution shows 
a more oscillatory surface pressure, especially on the lower surface, than other two solutions. The near-trailing-edge 
maximums of the lower-surface pressure in the CFL3D and TAU solutions are comparable and larger than that in the 
FUN3D solution. CFL3D and TAU show a small area of negative skin friction in the immediate vicinity of the trailing 
edge indicating flow separation; FUN3D shows no flow separation. Although not shown, FUN3D solutions on coarse 
grids also have some flow separation, but it goes away with grid refinement. 
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(a) Lift. 


(b) Pitching moment. 




(c) Viscous drag. 


(d) Pressure drag. 


Figure 9. Family II: Variation of forces and moment for the full airfoil. 
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(a) Lift. 


(b) Pitching moment. 




(c) Viscous drag. 


(d) Pressure drag. 


Figure 10. Family II: Variation of forces and moment constrained to the trailing edge region; 0.899166466843597 < x < 1.0. 
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(a) Lift. 


(b) Pitching moment. 




(c) Viscous drag. 


(d) Pressure drag. 


Figure 11. Family II: Variation of forces and moment constrained to the leading-edge region; 0.0 < x < 0.100177952877727. 
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Figure 12. Surface pressure and skin friction; moderate zoom; Family II grids. 
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(a) Cp at leading edge. 


(b) Cf jX at leading edge. 



(c) Cp at trailing edge. 



Figure 13. Surface pressure and skin friction; super zoom; Family II grids. 
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E. Grid Convergence at Different Locations 


This section provides a detailed description of the reference solutions computed on Family II grids. CFL3D, FUN3D, 
and TAU solutions are shown at several locations near the trailing edge and in the wake. Figures exhibit convergence 
and variations of the pressure, velocities and the turbulence variables. 

Figures 14-18 describe variations in the pressure coefficient near the trailing edge. Figure 14 shows the pressure 
variation along the vertical line corresponding to x = 1.001. This is a wake location in a close proximity to the 
trailing edge. The solutions computed by different codes and on different grids in the family are indistinguishable. 
This invariance indicates that an accurate pressure profile at this location can be computed on relatively coarse grids. 
Figures 15 and 16 present variations along the vertical line corresponding to x = 0.999; the variations corresponding 
to the upper (Fig. 15) and lower (Fig. 16) surfaces are shown separately. All solutions appear to be converging as 
grids are refined. CFL3D converges monotonically on the upper and lower surfaces. FUN3D and TAU converge 
monotonically on the upper surface. On the lower surface on the finest grid, the FUN3D solution appears to change 
the convergence direction. The coarse-grid TAU solution on the lower surface crosses the TAU solutions on finer 
grids; those finer-grid TAU solutions converge monotonically. On the upper surface, there is a noticeable difference 
between the CFL3D surface pressure and those of either FUN3D or TAU. The plots quickly become indistinguishable 
away from the surface. On the upper surface, FUN3D and TAU show larger variations in grid refinement than CFL3D 
solutions overall. However, the FUN3D and TAU solutions have a small variation on the 3 finest grids; the coarsest grid 
makes the solution variation large. On the lower surface, CFL3D shows larger variations between solutions computed 
on different grids than other two codes. Figures 17 and 18 show the pressure variation in the horizontal direction. 
Figure 17 shows variation behind the trailing edge (along the line z = 0); all plots are indistinguishable. Figure 18 
shows variations along the line 2 ? = 0.00008 located slightly above the trailing edge. Even though the finest grid 
solutions computed by all three codes are close to each other, all codes show significant variations between solutions 
computed on different grids. This variability indicates that finer grid resolution is required to accurately represent the 
local solution. 

Figures 19-21 show vertical variations of the horizontal velocity component, u, near the trailing edge. All plots 
are (almost) indistinguishable. Larger variations in u are shown in Fig. 22 in the wake region along the vertical line 
corresponding to x = 10. The variations are significant between solutions computed on different grids — the solutions 
on coarser grids do not resolve the wake velocity profile sufficiently. The variations between codes on the same grids 
are small. Figures 23 and 24 show variations of the u velocity in the horizontal direction near the trailing edge. The 
plots of the wake profile along z = 0 are indistinguishable. The variations along the 2: = 0.0008 line are significant 
between solutions computed by different codes and on different grids. The coarse-grid solutions indicate the presence 
of a reverse flow near the location corresponding to x ~ 0.9996. The reverse-flow pocket is larger in the CFL3D 
solution than in the FUN3D and TAU solutions. The reverse-flow pocket disappears on the finest grid for all solutions. 

Figures 25-30 demonstrate variation of the vertical velocity component, w, near the trailing edge and in the wake. 
Lines showing vertical variation near the trailing edge (Figs. 25-27) and the horizontal variation behind the trailing 
edge (Fig. 28) are indistinguishable. Large variations of ic-velocity are observed between solutions on different grids 
in the wake region along the line x = 10 (Fig. 29) and near the trailing edge along the line 2? = 0.00008 (Figs. 30). 
Variations between solutions computed with different codes on the same grids are small. 

Variations of the eddy viscosity near the trailing edge and in the wake are shown in Figs. 31-36. All eddy- 
viscosity plots in the near-trailing-edge region are close to each other. There are some visible differences at the edge 
of the boundary layer shown in Figs. 31-33. CFL3D shows a larger variation between solutions on different grids than 
FUN3D. The TAU solution on each grid shows a small local oscillation of the eddy viscosity at 2 ? > 0.047 (Figs. 31 
(c) and 32 (c)) and another one at 2: < —0.017 (Fig. 33 (c)); other codes show a smooth transition to zero in these 
regions. The explanation for this oscillation is that TAU uses a central difference scheme with a small matrix-valued 
artificial dissipation for inviscid fluxes, while the other two codes use upwind-biased schemes for inviscid fluxes. Note 
that the amplitude and footprint of the oscillation quickly decrease with grid refinement. 

Similar to the meanflow characteristics, eddy viscosity has a significant variation in the wake (Fig. 34) and a small 
variation behind the trailing edge (Fig. 35). Note that CFL3D solutions on the two finest grids are not shown in Fig. 35 
(a). With the SA model residuals at the level of 10“ 7 , the CFL3D eddy- viscosity profiles along the cut 2? = 0 were 
still changing, albeit very slowly. The computations on the two finest grids were stopped before converged eddy- 
viscosity profiles have been achieved. The observed wake variations are mainly between eddy viscosity computed on 
different grids. Similar profiles are obtained by different codes on the same grids. In distinction from the meanflow 
characteristics, the grid variations of eddy viscosity near the trailing edge along the line 2: = 0.00008 (Fig. 36) are 
small and plots are very similar for all solutions. Although not shown, eddy- viscosity convergence is very sensitive to 
the approximation order for the convection term in the SA model equation. CFL3D solutions computed with the first- 
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order approximation showed a significant deterioration of accuracy and convergence for the eddy- viscosity profiles 
near the trailing edge and in the wake. 

Figure 37 shows variations of the non-dimensional Spalart turbulence variable in FUN3D solutions. Figures 37 (a) 
and (b) zoom to two locations, the boundary-layer edge and the wake edge, where negative values of the turbulence 
variable are observed. Near the boundary-layer edge, the coarser-grid solutions show negative turbulence values of 
larger amplitude than fine-grid solutions. On finer grids, the area with negative turbulence variables is significantly 
reduced. However, the number of nodes with negative turbulence values appears approximately constant on all grids. 
Near the wake edge, the area of negative turbulence variables decreases on finer grids, but the number of affected 
nodes does not decrease. The amplitude of the negative turbulence initially increases in grid refinement, but decreases 
on the finest grid. 



(a) CFL3D. 



(b) FUN3D. 



0,04 


(c) TAU. 


Figure 14. Cp variation behind the trailing edge along the line x = 1.001. 





Figure 15. Cp variation along the line x = 0.999 over the upper surface. 


17 of 50 


American Institute of Aeronautics and Astronautics 



(a) CFL3D. 




(c) TAU. 


Figure 16. Cp variation along the line x = 0.999 under the lower surface. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 17. Cp variation behind the trailing edge along the line z = 0. 





(a) CFL3D. (b) FUN3D. (c) TAU. 

Figure 18. Cp variation near the trailing edge along the line z = 0.00008. 
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(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 19. u-velocity variation behind the trailing edge along the line x = 1.001. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 20. u-velocity variation along the line x = 0.999 over the upper surface. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 21. u - velocity variation along the line x = 0.999 under the lower surface. 
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u/a^j u/a rflP u/a^, 


(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 22. u- velocity variation in the wake along the line x = 10. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 23. u-velocity variation behind the trailing edge along the line z = 0. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 24. u-velocity variation near the trailing edge along the line z = 0.00008. 
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(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 25. ^-velocity variation behind the trailing edge along the line x = 1.001. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 26. lovelocity variation along the line x = 0.999 over the upper surface. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 27. lovelocity variation along the line x = 0.999 under the lower surface. 
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(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 28. lovelocity variation behind the trailing edge along the line z = 0. 



(a) CFL3D. (b) FUN3D. (c) TAU. 


Figure 29. velocity variation in the wake along the line x = 10. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 30. ^-velocity variation near the trailing edge along the line z = 0.00008. 
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(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 31. Eddy- viscosity variation behind the trailing edge along the line x = 1.001. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 32. Eddy- viscosity variation along the line x = 0.999 over the upper surface. 
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(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 33. Eddy- viscosity variation along the line x = 0.999 under the lower surface. 



-o.oi 

o 

h) 

0 02 


*003 


0 50 100 150 200 250 

JVM™ 



100 150 

JVM™ 


- 0.01 


23 of 50 


American Institute of Aeronautics and Astronautics 




u 

a 


2 
1.B 
1.6 
1.4 
1.2 
1 

0 TOOO 2005 



(a) CFL3D. 




(b) FUN3D. (c) TAU. 


Figure 34. Eddy- viscosity variation in the wake along the line x = 10. 





(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 35. Eddy-viscosity variation behind the trailing edge along the line z = 0. 



(a) CFL3D. 


(b) FUN3D. 


(c) TAU. 


Figure 36. Eddy- viscosity variation near the trailing edge along the line z = 0.00008. 


24 of 50 


American Institute of Aeronautics and Astronautics 






0,08 
0.075 
0,07 
0.065 
^ 0.06 
0.055 
0,05 
0.045 
0,04 

-5 0 5 10 15 20 

SA-neg turbulence variable {nondim) 











G finest grid 

- -A t level coarser 

v 2 levels coarser 

p 3 levels coarser 
















1 

/ 

r 





E 

1 

b* 
















1 1 



SA-neg turbulence variable (nondim) 


(a) x = 0.999, upper surface, zoom to boundary layer edge. 


(b) x = 10, zoom to wake edge. 


Figure 37. FUN3D: Variation of the non-dimensional Spalart turbulence variable. 


F. Effects of Discretization Methods and Grid Elements 

Effects of variations in the FUN3D discretization methods and grid elements are studied for grids of Family I and 
Family II. The discretization methods vary in the approximation order of the convection term in the turbulence-model 
equation and the reconstruction method for the inviscid terms in the meanflow equations. 

The baseline results correspond to a nominally second-order approximation of meanflow and turbulence-model 
equations. The meanflow variables are reconstructed from a node to the edge midpoints using the k = 0.5 MUSCF 
scheme. In the direction approximately tangential to the airfoil surface and horizontal in the wake, least-squares (lsq) 
gradients are used for reconstruction. The lsq gradients at a node are constructed from an unweighted least- squares 
linear fit to the primitive variables at the neighboring nodes. In the direction approximately normal to the airfoil 
surface and vertical in the wake, the gradients are computed using an implicit mapping along the structured grid lines 
approximately following this direction. Such gradients are denoted as mapped-n gradients. The second-order accurate 
discretization of the SA model convection term always uses lsq gradients for reconstruction of the Spalart turbulence 
variable from the node to the edge midpoint. The reconstruction corresponds to the k = 0 scheme. The baseline 
results on quadrilateral grids of Family II are expected to be the most accurate of the results presented. 

Figures 38-40 compare convergence of the lift, pitching moment, viscous drag, and pressure drag coefficients. 
Each plot is characterized by four parameters: grid Family (I or II), grid element type (Q for quadrilateral elements or T 
for triangular elements), convection approximation order in the SA equation (1st or 2nd), and the type of gradients used 
for reconstruction in the normal direction (mapped-n or lsq). For example, the baseline discretization plot is designated 
as I: Q: 2nd: mapped-n. Figure 38 shows convergence plots for different discretization methods on quadrilateral grids. 
The grid element designation (Q) is omitted as it would be the same for all plots on the figure. The lift and pitching 
moment computed on grids of the same family with different discretization methods differ little in comparison to 
the differences between solutions on grids of different families. Note that errors introduced by poor trailing-edge 
resolution (on Family I grids) negate effects of more accurate discretization methods. In particular, lift and pitching 
moment coefficients computed on Family I grids with mapped-n gradient reconstruction appear less accurate than the 
coefficients computed with lsq gradient reconstruction. The viscous drag coefficient computed from solutions with the 
first-order approximation in the SA equation is significantly lower than the coefficient computed from a solution with 
a second-order approximation, independent of the grid family. In pressure drag, the variations due to decreasing the 
accuracy of the turbulence-model convection term and the meanflow reconstruction method are in opposite directions. 
Thus, the pressure drag coefficients computed from the presumably most and least accurate solutions are similar on 
the third-finest grid. In either viscous or pressure drag, the total variation across all solutions on the finest grid is 
approximately 0.1 count. 

Convergence of the aerodynamic forces and moment on triangular grids of Family I and Family II is shown in 
Fig. 39 and Fig. 40, respectively. The results on quadrilateral grids of Family I and II are shown for reference. 
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Figure 38. Effect of variations in discretization method for quadrilateral grids. 
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Approximation for the convection term in the S A model equation has a significant effect on solution accuracy on grids 
with triangular elements. The lift computed with the first-order approximation is reduced and the pitching moment is 
increased in comparison with the corresponding quantities computed using a second-order approximation. Solutions 
on triangular grids of Family I produce higher lift and pressure drag and a lower pitching moment than solutions on 
corresponding quadrilateral grids. Similar to computations on quadrilateral grids, the order of the convection-term 
approximation in the SA equation is the major factor affecting accuracy of the viscous drag coefficient. The viscous 
drag is significantly lower with the first-order approximation than with a second-order approximation independent of 
the grid elements, family, and meanflow flux reconstruction scheme. The sensitivities to triangular elements are more 
pronounced for Family I grids than for Family II grids. The solutions on triangular Family II grids computed with a 
second-order accurate SA-model convection term are very similar to the baseline solutions on quadrilateral Family II 
grids. 
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Figure 39. Triangular grids of Family I: effect of discretization methods. 
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Figure 40. Triangular grids of Family II: effect of discretization methods. 


28 of 50 


American Institute of Aeronautics and Astronautics 


IV. Flat Plate Configuration 


A grid convergence study for a turbulent flow through a cascade of finite flat plates is presented. FUN3D is used to 
establish an accurate reference solution and to assess effects of streamwise grid resolution near geometric singularities 
(i.e., the leading and trailing edges) on convergence of turbulent-flow solutions. 

A. Test Case Description 



Figure 41. Boundary conditions for the finite flat plate geometry; 81 x 25 grid is shown. 


This test case corresponds to the “2D Finite Flat Plate” case defined under the “Cases and Grids for Turbulence 
Model Numerical Analysis” section at the TMR website. The coordinate system is defined using non-dimensional 
units; x and z are streamwise and normal (vertical) coordinate directions, respectively. The boundary conditions 
are shown schematically in Fig. 41. The flat plate is located at the bottom of the domain (0 < x < 2,z = 0), 
similarly to the “2D Zero Pressure Gradient Flat Plate” case defined under the “Turbulence Model Verification Cases 
and Grids” section at the TMR website. For the current study, the computational domain (— 2 < x < 4, 0 < 
z < 4) has been extended upward, upstream, and downstream of the plate, creating a trailing edge inside of the 
computational domain at xte = 2, zte = 0. The specific placement of the top, upstream, and downstream boundaries 
has been chosen to ensure that the drag coefficient is within 0.02 counts of the coefficient computed on a domain 
with such boundaries placed at infinity (see Sec. IV.F below). The top boundary condition has been changed to a 
symmetry condition to avoid specifying the external state along the top of the computational domain. Thus, the test 
case corresponds to a cascade of finite flat plates separated by a distance of eight in the vertical ^-direction. Adiabatic 
no-slip boundary conditions are applied at the plate surface (0 < x < 2, z = 0). Constant total pressure boundary 
conditions, corresponding to T t /T re f = 1 + 0.2M^ and P t /P re f = (T t /T re j) 3,5 , are applied at the upstream 
boundary (x = —2, 0 < z < 4). Constant pressure boundary conditions, corresponding to P/P re f = 1, are applied at 
the downstream boundary {x = 4, 0 < z < 4). Symmetry conditions are applied at the top (— 2 < x < 4, z = 4) and 
bottom (—2 < x < 0 and 2 < x < 4, z = 0) boundaries. 

FUN3D computations have been performed on a series of nested, stretched, rectilinear grids ranging from the 
2, 561 x 769 (finest) grid to the 21 x 7 (coarsest) grid. The numbers in the grid notation indicate the points in the 
streamwise and normal directions, respectively. Each coarser grid is exactly every-other-line of the next finer grid. 
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The grids are stretched in the wall-normal direction and are clustered near the leading and trailing edges. The normal 
spacing of the finest grid at the wall is 2.5 x 10 -7 , corresponding to a non-dimensional boundary-layer spacing of 
z + =0.1 at the middle of the plate (x « 1). The x-directional clustering near the leading edge is set by specifying 
a local aspect ratio ( ARle ) of the grid. The grids are symmetric about x = 1 and thereby the trailing-edge aspect 
ratio is the same as at the leading edge. The recommended value, ARle = 1, is used by default, although some 
computations on grids with ARle = 1200, which is more typical of current practice, are reported in Sections IV.B 
and IV.C below. Figure 41 shows a view of the 81 x 25 grid. 

The flow conditions correspond to = 0.2 and Re = 5M based on a unit length of the grid. The body 
reference length is two units. Thus, Re = 5M at the middle of the plate at x = 1 and Re = 10M at the trailing 
edge of the plate at x = 2. The SA-neg version of the Spalart-Allmaras turbulence model variable is used, although 
at convergence there are no negative values of the turbulence variable. The farfield value of the Spalart turbulence 
variable is v far field = 3^. The Prandtl number is taken to be constant at Pr = 0.72, and the turbulent Prandtl 
number is taken to be constant at Pr t = 0.9. The molecular viscosity is computed using Sutherland’s Law. 

B. Drag and Maximum Eddy Viscosity 




(a) Drag coefficient. 


(b) Maximum eddy viscosity. 


Figure 42. Grid convergence of drag and maximum eddy viscosity. 


The convergence plots of the drag coefficient and the maximum eddy viscosity are shown in Fig. 42 for ARle = 1 
and ARle = 1200 versus h = yJ\/N. The drag scale is quite fine, spanning only 0.1 drag count. Both the 
drag coefficient and the maximum eddy viscosity show less variation in grid convergence with ARle = 1 than 
with ARle = 1200. The maximum eddy viscosity converges linearly in h for each aspect ratio, indicating first- 
order convergence. Close examination shows that the maximum eddy viscosity occurs just above the wake centerline 
(see Sec. IV.D below). The first-order variation of maximum eddy viscosity is believed to be a boundary effect, as 
grid convergence at a fixed location upstream of the trailing edge is approximately second order. In contrast, drag 
convergence is first order on the three finest grids for ARle = 1 and less than first order for ARle = 1200. From 
the boundary layer theory and the numerical results here, the skin friction in the leading edge varies as 0( 1/y/x). 
The drag integration routine is a trapezoidal second-order integration, so drag convergence is expected to be first order 
even if the skin friction values were exact. 

To investigate convergence of local drag contributions, three sections on the surface are chosen: near the leading 
edge, in the middle of the plate, and near the trailing edge; the sections are defined in Table 3. The separator nodes, 
(x,z) = (0.107267441655523,0) and (x,z) = (1.89273255834448,0), are present on the finest four AR L e = 1 
grids used in the investigation. Convergence of drag within each of the sections is shown versus h p in Fig. 43 for 
various choices of p. The dashed lines in the figure are linear fits for the finest two grids. In the leading-edge section, 
the drag convergence order is slightly less than first, as would be expected. In the middle-plate section, the drag 
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(a) Leading-edge section. 


(b) Middle section. 


(c) Trailing-edge section. 


Figure 43. Grid convergence of drag within different sections. 


convergence order is high, close to p = 3.3, which is unexpected. In the trailing-edge section, the drag convergence 
order is slightly greater than first. 

Table 3. Sections of the plate. 

Leading Edge Middle Plate Trailing Edge 

0 < x < 0.107267441655523 0.107267441655523 < x < 1.89273255834448 1.89273255834448 < x < 2 


C. Skin Friction, Surface Pressure, and Boundary-Layer Profile 

The skin friction variation near the leading edge is shown in a global view (Fig. 44 (a)) and in a local view near 
x = 0.000101 (Fig. 44 (b)). Four finest grids with ARle = 1 are used: “grid 1” designates the finest 2, 561 x 769 
grid, “grid 4” designates the 321 x 97 grid. An approximate analytic fit to the skin friction variation corresponding to 
0. 0002503a; -0 - 52 is shown, with the constant selected to match the infinite-grid extrapolated value at x = 0.000101 
assuming second-order grid convergence. In Fig. 44 (a), the maximum discretization error actually grows on finer 
grids. The maximum relative discretization error, defined as the local discretization error divided by the local solution 
value, reduces quite slowly in grid refinement. In Fig. 44 (b), the skin friction plots on the two finest grids are 
indistinguishable. The skin friction at x = 0.000101 is shown in Fig. 45 versus h and h p to assess grid convergence. 
The dashed lines in the figures fit the results on the two finest grids. Grid convergence exhibits an order very close to 
1.6 on the three finest grids. 

Skin-friction convergence in the middle of the plate at x = 1.0 is shown in Fig. 46 for ARle = 1 and ARle = 
1200. The ARle = 1 results show an apparent convergence order of 2.5; the ARle = 1200 results show less-than- 
first convergence order. Somewhat surprisingly, the ARle = 1200 results have smaller discretization error levels on 
coarser grids; the errors are extremely small for all solutions. 

The skin friction variation near the trailing-edge region is shown in a global view (Fig. 47 (a)) and near x = 1.95 
(Fig. 47 (b)). Results with both ARle = 1 and ARle = 1200 are shown. The authors could not find a simple function 
to characterize the variation of the skin friction in this range; the function 0.0027001033x - ° 14 + 0.000007/(2 — x) 0,5 
is shown to indicate that the skin friction is singular near the trailing edge. Assuming that the ARle = 1 solution 
on grid 1 is reasonably close to the exact solution, large errors in ARle = 1200 solutions are evident over the four 
grid points closest to the trailing edge and the maximum ARle = 1200 discretization error decays very slowly, if at 
all. The ARle = 1 solutions also exhibit the largest errors at the four grid points closest to the trailing edge. Grid 
convergence at the trailing edge is not analyzed in detail, but can be inferred from Fig. 47 (b). The ARle = 1 results 
exhibit low discretization errors and a better-than- first grid convergence order, while the ARle = 1200 results show 
larger discretization errors and an apparent less-than-first grid convergence order. 

Convergence of the pressure coefficient, C p , near the leading and trailing edges is shown in Fig. 48. The pressure 
appears discontinuous at both the edges along the line tangential to the plate. Grid convergence is slow at the grid 
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(a) Global view. 


(b) Near x = 0.0001. 


Figure 44. Skin friction convergence in leading-edge region; ARle = 1. 



Figure 45. Skin friction grid convergence at x = 0.000101; ARle = 1 
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Figure 46. Skin friction grid convergence at £ = 1.0. 


tfj gndl 




(a) Global view; cl = 0.0.0027001033; c2 = 0.000007. 


(b) Near x = 1.95. 


Figure 47. Skin friction convergence near trailing edge. 
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points immediately adjacent to the edges but, at a fixed distance away from the edges, convergence is approximately 
second order. Contours of the pressure coefficient near the leading and trailing edges are shown in Figs. 49 (a) and (b), 
respectively. Near both the edges, the pressure coefficient varies smoothly along rays emanating from the edges. 




(a) Near leading edge. 


(b) Near trailing edge. 


Figure 48. Surface pressure coefficient. 



cp 

0 4 

0.3002 

0.3001 

03106 

0.3208 

0301 

02012 

02014 

024 16 

02218 

0202 

01022 

01001 

01426 

01228 

0103 

00032 

0.0631 

00136 

00238 

0001 


X 



-5E-05 - 


1.99995 


2 00005 2 0001 


0.00015 


0.0001 


5E-05 

N 


- 0.02 
-0.021 
' 0.022 
0023 
0.021 
-0.02b 
■0026 
- 002 f 
0.026 
0.029 
-OUB 


(a) Near leading edge. 


(b) Near trailing edge. 


Figure 49. Contours of pressure coefficient. 

Figure 50 shows the boundary layer profiles of the horizontal velocity, u , in the middle of the plate at x = 1. In 
the global view with a logarithmic z-scale, the velocity profiles computed on different grids are indistinguishable from 
each other. In a detailed view near 2 « 0.001, grid convergence is qualitatively second order. Similarly, in the global 
view (Fig. 51 (a)), the eddy viscosity profiles are indistinguishable from each other; and the detailed view near the 
edge of the boundary layer (Fig. 51 (b)) reveals grid convergence that is qualitatively better than first order. 
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(a) Global view. 


(b) Near middle of logarithmic region. 


Figure 50. Velocity profile at x = 1. 



(a) Global view. 


(b) Near edge of boundary layer. 


Figure 51. Eddy viscosity profile at x = 1. 
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D. Downstream of Trailing Edge 


Figures 52 and 53 show grid convergence of the horizontal velocity, u, downstream of the trailing edge. Away from 
the trailing edge, the horizontal velocity increases as (x — xte ) 0 ' 08 and the velocity distributions on different grids 
are nearly indistinguishable from each other in the global view (Fig. 52 (a)). Near the trailing edge (Fig. 52 (b)), the 
maximum errors are reducing slowly as the grid is refined. At a small distance from the trailing edge, the velocity 
variation slope is changed to (x — xte) 0A 3 > and the errors converge with the rate between first and second order. 
Figures 53 and 54 illustrate convergence at locations corresponding to close proximity to the trailing edge (( x—xte ) ~ 
10 -4 ) and at a location further in the wake (x « 3). At both locations the apparent convergence order is 1.6. Note 
that the turbulent-flow velocity profile near the trailing edge is similar to the profile of a reentrant-corner solution for 
a pure-diffusion equation described in the Appendix, but the observed convergence rate is higher. The pure-diffusion 
solution exhibits a square-root behavior near the singularity and first-order convergence at any fixed interior location. 




Figure 52. Velocity downstream of the trailing edge. 


Figures 55 and 56 show grid convergence of the pressure (C p ) and eddy viscosity (fi t ) downstream of the trailing 
edge. The pressure distributions (Fig. 55 (a)) on different grids are nearly indistinguishable from each other, except in 
an immediate vicinity of the outflow boundary, where the specified-pressure boundary condition forces sharp solution 
variations over a few grid points adjacent to the boundary. Fig. 55 (b) indicates that, in the trailing-edge vicinity, 
maximum errors in pressure are at the grid points nearest the edge. Note grid convergence in Fig. 55 (b) is the same 
as in Fig. 48 (b), but shown with a logarithmic scale for the abscissa. At a fixed distance away from the trailing edge, 
grid convergence is qualitatively second order. 

The eddy viscosity distributions behind the trailing edge computed on different grids are nearly indistinguishable 
from each other (Fig. 56 (a)). Although not shown, examination of the eddy viscosity at a fixed distance from the 
trailing edge shows better- than-first-order grid convergence. 

Figures 57 and 58 show the wake profiles of the horizontal velocity and eddy viscosity at x = 3. In the global 
view, the velocity profiles computed on different grids are indistinguishable from each other except near the edge of 
the wake. In a detailed view near the edge of the wake, grid convergence is qualitatively second order. Similarly, the 
eddy viscosity profiles from the three finest grids are close to each other in the global view with small deviations near 
the edge of the wake; and the detailed view near the maximum-eddy- viscosity location reveals grid convergence that 
is qualitatively second order. 

E. Iterative Convergence 

Figures 59 (a) and (b) illustrate iterative convergence of the Full Multigrid (FMG) solver 35 on grids with ARle = 
1200 and ARle = 1, respectively. Convergence of the L 2 norms of the meanflow and turbulence-model residuals 
and the drag coefficient is shown. Only four grids are used in the FMG process. A Full Approximation Scheme (FAS) 
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(a) Local view. 


(b) Convergence order. 


Figure 53. Velocity convergence at x — xte ~ 10 4 . 




(a) Local view. 


(b) Convergence order. 


Figure 54. Velocity convergence at x « 3. 
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(a) Global view. (b) Near trailing edge. 


Figure 55. Pressure downstream of the trailing edge. 



(a) Global view. 


(b) Near trailing edge. 


Figure 56. Eddy viscosity downstream of the trailing edge. 
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(a) Global view. 


(b) Near edge of wake. 


Figure 57. Wake velocity profile at x = 3. 



(a) Global view. 


(b) Near maximum eddy viscosity. 


Figure 58. Wake eddy viscosity profile at x = 3. 
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nonlinear V-cycle, 35 FAS (2, 2), with 2 relaxations preceding and 2 relaxations following the coarse-grid correction, 
is used. Each relaxation uses a tightly-coupled formulation and one multicolor point-implicit sweep and one line- 
implicit sweep over the entire domain at each multigrid level. The coarsest grid (level 4) takes just over 100 cycles 
(400 relaxations) to converge residuals to the residual tolerance of 10 -12 . The finer grids take about 20 FAS cycles to 
converge residuals to the tolerance. The drag coefficient converges faster on grids with ARle = 1 than on grids with 
ARle = 1200; specifically, the level of discretization error is reached in just a few cycles on grids with ARle = 1 
and within ten cycles on grids with ARle = 1200. This faster convergence of the drag coefficient is attributed to a 
better initial approximation provided by the FMG solver from a coarser grid with ARle = 1. 




(a) AR le = 1200. (b) AR LE = 1. 

Figure 59. Iterative convergence of residuals and drag coefficient versus multigrid FMG cycle for 4 finest grids. 


F. Variation of the Farfield Boundary Locations 

The effects of upper and upstream/downstream boundary locations have been studied parametrically and the results 
are shown in Fig. 60. For these studies, the aspect ratio of the grids was ARle = 1200. The grids are symmetric about 
x = 1, so the distances from the plate to the upstream and downstream boundaries are the same. The drag coefficient 
and maximum eddy viscosity on the finest grid are shown for varying locations. In Fig. 60 (a), the distance from the 
plate to the upstream and downstream boundaries is fixed at 1 (— 1 < x < 3), and the upper boundary is shifted. Both 
the drag and maximum eddy viscosity vary linearly with respect to the inverse of distance from the plate to the upper 
boundary. The drag coefficients computed on the domains with upper boundary located at z = 4 and 2 = 1 are 0.015 
and 0.055 counts lower, respectively, than the drag coefficient extrapolated to the limit of the boundary at z = 00 . 
Corresponding changes to the maximum eddy viscosity are less than 1%. 

In Fig. 60 (b), the upper boundary fixed at z = 4, and the inflow/outflow boundary locations are shifted. The drag 
varies linearly with respect to the inverse of distance from the plate to the downstream (or upstream) boundary location. 
The variation in the computed drag coefficient is smaller (< 0.01 drag count) than the variation due to changes 
in the upper boundary location. The maximum eddy viscosity increases considerably as the upstream/downstream 
boundaries are moved farther from the plate, slightly faster than the inverse of the distance to the boundaries. Based 
of these results, the upper boundary location was chosen as z = 4 and the horizontal extent of the domain was chosen 
as —2 < x < 4 to provide the distance of 2 from the plate to the downstream (and upstream) boundary. 


V. Concluding Remarks 

A detailed grid convergence study has been conducted to establish accurate reference solutions corresponding to a 
one-equation linear eddy-viscosity Spalart-Allmaras (SA) turbulence model for two-dimensional (2D) turbulent flows 
around the NACA 0012 airfoil and a flat plate configuration. The investigation of the NACA 0012 airfoil involved 
three widely used computational fluid dynamics (CFD) codes, FUN3D (NASA), CFL3D (NASA), and TAU (DLR), 
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(a) Variation in upper boundary location; distance from plate 
to downstream boundary of 1 . 



V V.c. U.-H- U.u U.o I 

1 /(Distance to Inllow/Outflow Boundary From Plate) 


(b) Variation in downstream boundary location; distance 
from plate to upper boundary of 4. 


Figure 60. Effect of varying boundary locations; AR LE = 1200; FUN3D finest grid values. 


and three families of uniformly refined structured grids with different density distribution. The following observations 
have been reported. 

1. Solutions computed on different grid families appear to converge to the same continuous limit, but exhibit 
different convergence characteristics. 

2. The grid resolution in the vicinity of geometric singularities, such as a sharp trailing edge, is the major factor 
affecting accuracy of discrete solutions, more prominent than differences in discretization schemes and/or grid 
elements. 

3. On grids from the family with an improved trailing-edge resolution, the solutions obtained with different codes 
are similar. Plotted on a global scale, the solutions on the finest grid are almost indistinguishable. Differences in 
the pressure and skin friction coefficients appear only in a narrow range within 0.001c distance from the trailing 
edge and near the minimum pressure location in the leading edge. Off-body profiles differ mostly in this narrow 
range of the trailing edge. 

4. The aerodynamic coefficients predicted by the three codes on the finest grid with 15M degrees of freedom and 
with an improved trailing-edge resolution show an impressive agreement. The code-to-code variations in the 
total drag are less than 0.1% (0.1 count), in the lift are less than 0.02%, and in the pitching moment are less than 
1 %. 

5. Even on such fine grids producing such accurate solutions, the asymptotic convergence order has not been 
established. 

Similar observations have been made for FUN3D computations performed for a turbulent flow around a cascade of 
flat plates: the solution accuracy and convergence have greatly benefited from an improved stream wise resolution near 
geometric singularities, i.e., the leading and trailing edges. Namely, solutions computed on a sequence of uniformly 
refined grids with a high resolution near the leading and trailing edges (i.e., with the local aspect ratio of ARle = 1) 
show significantly less variation in grid refinement than solutions computed on grids with the same number of degrees 
of freedom but with a lower resolution near the edges (the local aspect ratio of ARle = 1200). On grids with 
ARle = 1 , the drag contributions from the leading and trailing-edge sections converge with first order. The drag 
contribution from the middle-plate section converges with an apparent order of 3.3. Skin friction convergence exhibits 
an order of 1.6 in a fixed location next to the leading edge, an order of 2.5 in the middle of the plate, and a better- 
than-first order near the trailing edge. The flow velocity at fixed locations behind the trailing edge converges with an 
apparent order of 1.6. The local eddy viscosity converges with an apparent second order in the middle of the plate, 
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near the trailing edge, and in the wake, but the maximum eddy viscosity converges with approximately first order. 
On grids with ARle = 1200, the skin friction errors are less than those with ARle = 1 on coarser grids but the 
grid convergence order is less than first order. On all grids, the maximum errors in skin friction and surface pressure 
coefficients always occur over a fixed number of grid points near the leading and/or trailing edges and decay slowly 
with grid refinement. 

Note that on the grids used in the study, convergence of turbulent-flow solutions near geometric singularities 
is significantly different from the solution for an elliptic (pure diffusion) equation. The elliptic-equation solution 
exhibits first-order convergence on uniformly refined grids at any fixed interior location. Observed convergence of the 
turbulent-flow solutions does not degrade to the same degree. This difference in convergence may indicate that, in 
spite of a high solution accuracy obtained on the grids used in the study, much finer grids are needed to observe the 
expected first-order asymptotic convergence order. This observation is discouraging to the prospects of realizing high 
asymptotic convergence orders on uniformly refined grids. The topic deserves further investigation. A corner (non- 
uniform) refinement strategy is shown in the appendix to recover design order convergence for the elliptic-equation 
solution. Its success in application to the pure-diffusion equation provides an indication that improved resolution near 
geometric singularities is essential to improve the accuracy of turbulent-flow simulations. 
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Appendix: Elliptic Equation with Reentrant Corner Singularity 

The long-known accuracy degeneration in solutions of elliptic equations on domains with reentrant comer singu- 
larities 11 is revisited. The reentrant corner geometry is represented parametrically by an exterior angle, uj, that can 
be selected to match the local trailing-edge geometry of the NACA 0012 airfoil. The computations shown in this 
appendix are for uj = 27 r, which corresponds to a cusped trailing-edge geometry and the most severe accuracy degra- 
dation. Elliptic equations describe diffusion phenomena and thus applies directly to Stokes flows. The relevance of 
the accuracy degeneration to high Reynolds number flows has not been studied. 

The exact solution, q , of the 2D Laplace equation, 


&xxQ H" dyyQ — 0? 

satisfying homogeneous Dirichlet boundary conditions for a reentrant corner situated on the left of the computational 
domain (as in Fig. 61 (a)) is given by the following, 

q = r a sin (aO L ), 

where the angle 6 L increases clockwise with a branch cut along 6 = 7r, 

e L = 7T - e. 

The position of the reentrant comer tip (also referred as the trailing edge by analogy with the NACA 0012 geometry) 
is (xte, zte ) = (1,0), the inverse relative exterior angle is a = tt/uj, and 


X = X — XtEj 
z = z - zte , 

6 = arctan2(z, x), 
r = y/ x 2 + z 2 . 

Dirichlet conditions, q = 0, are applied for x < xte along the radial lines 6 = it — (2i r — uS)/2 and 0 = —7 r + (27 r — 
cc)/2. The exact solution contours with uo = 2e are shown in Fig. 61 (a). For completeness, the exact solution with the 
reentrant corner situated on the right of the domain, corresponding to the geometry most often cited in the literature, 35 
is given by the following, 


q = r a sin (a0 R ), 
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where the angle 6 R increases counter-clockwise with a branch cut along 6 = 0, 


e R 


6 + 2tt if z < 0, 

6 otherwise. 


The theoretical estimates 11 for convergence orders of discretization error, e, on uniformly refined Cartesian grids 
are 


Moo*h & , 

Here, 1 1 • 1 1^ and 1 1 • | |i are the and L\ norms, respectively. 

The lowest convergence order is predicted for a cusped trailing edge, which is a zero-angle reentrant corner (u = 27 r, 
a = 0.5), 


||e||oo « Vh, m 

Iklli ~ h. 

This is a significant global degradation of second-order convergence expected for nominally second-order discretiza- 
tion schemes and smooth solutions. 

To demonstrate and repair this degradation, several sequences of rectilinear grids are generated on the domain 
(x,z) G [0,2] x [—1,1] centered about (xte, zte)- Solutions are computed on a sequence of uniformly-refined 
isotropic Cartesian grids as well as on sequences of non-uniformly refined rectilinear grids with additional degrees of 
freedom added near the reentrant corner. The latter refinement is reference as corner refinement to distinguish it from 
uniform refinement. 

A 2D rectilinear grid is derived from a one-dimensional (ID) primal mesh generated on the interval 0 < x < 1. 
This ID mesh is mapped symmetrically on the interval 1 < x < 2 and also on the z-axis. The 2D mesh is constructed 
as the tensor product of the ID x- and z-directional meshes (see Fig. 61 (b) and (c)). A sequence of uniform nested 
ID meshes on the primal interval leads to a sequence of uniformly refined isotropic Cartesian grids. 





(a) Exact solution. The reentrant-corner 
boundary is the black line on the left side 

(x < 1, z = 0). 


(b) Corner-Refinement grid (level 4; 2 
seed cells). 


(c) Corner-Refinement grid near reen- 
trant corner (level 4; 2 seed cells). 


Figure 61. Reentrant corner: exact solution and grids. 


The starting coarse grid for a sequence of corner-refined rectilinear grids is one of the Cartesian grids characterized 
by the number mesh spacings (cells) over the interval 0 < x < 1; this number is referred as the number of seed cells. 
One step of the grid refinement is defined as follows. First, all coarse cells except the reentrant-comer cell (the cell 
attached to the (xte, z te ) node) are divided into two equal fine cells. The coarse reentrant-corner cell undergoes the 
corner refinement. It is subdivided 4 times; in each subdivision, only the local cell closest to the reentrant corner is 
divided into two equal cells. Thus, the original coarse reentrant-corner cell has been divided into four cells, while each 
of the other coarse cells has been divided into two cells. If the size of the coarse reentrant-corner cell is defined by 
5 = 1, then the sizes of the four fine cells are s = {s/2, s/4, s/8, s/ 16, s/16}, where the two smallest spacings are 
nearest to the reentrant corner. Note that the mesh size away from the reentrant-corner has been uniformly reduced by 
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factor 2; the mesh size next to the reentrant corner has been reduced by factor 16. The further mesh refinement on the 
interval 0 < x < 1 is done recursively. 

The mesh size, h, in the theoretical estimate of Eq. 1 is the local mesh size in a vicinity of the reentrant comer. 
In the comer refinement, the square root of the local mesh size is proportional to the square of the mesh size in the 
uniformly refined part of the domain. This corner refinement is expected to recover second-order convergence of 
discretization errors for the worst case of uo = 2i r. 

The number of seed cells effectively determines the physical extent of the grid affected by the corner refinement; 
the rest of the computational domain undergoes the uniform grid refinement. A grid started from 2 seed cells (level 1) 
and refined through 3 steps of comer refinement (to level 4) is shown in Fig. 61 (b) and (c). The grid differs from a 
corresponding uniformly refined grid over « 44 percent of the domain. For example, the upper left quadrant defined 
by x < 0.75 and z > 0.25 is uniformly refined. A corner-refinement grid started from 16 seed cells (not shown) differs 
from a uniformly refined over a much smaller region (« 6 percent) of the physical domain; the upper left quadrant 
defined by x < 0.9675 and 2 > 0.0325 is uniformly refined. 

The degrees of freedom on the uniform and comer-refinement grids are shown in Tables 4 and 5 for grid sequences 
started from 2 and 16 seed cells, respectively. The degrees-of-freedom increase from level to level is also tabulated for 
the comer refinement. Although the increase factor is quite large for coarse-grid levels, the factor asymptotes to the 
factor of four, which is typical for the uniform refinement. Note that the uniform-refinement grid levels in the tables 
are not synchronized; level 2 from Table 4 corresponds to level 1 from Table 5. 

Table 4. Number of degrees of freedom on the uniform and corner-refinement grids with 2 seed cells. Number in parentheses is the 
degrees-of-freedom increase factor from the previous level. 


Grid Fevel 

Uniform-Refinement 

Corner-Refinement 

8 (coarsest) 

27 

27 

7 

85 

232 (8.59) 

6 

297 

1242 (5.35) 

5 

1105 

5662 (4.56) 

4 

4257 

24102 (4.26) 

3 

16705 

99382 (4.12) 

2 

66177 

403542 (4.06) 

1 (finest) 

263425 

1626262 (4.03) 


Table 5. Number of degrees of freedom on the uniform and corner-refinement grids with 16 seed cells. Number in parentheses is the 
degrees-of-freedom increase factor from the previous level. 


Grid Fevel 

Uniform-Refinement 

Corner-Refinement 

4 (coarsest) 

1105 

1105 (NA) 

3 

4257 

5076 (4.59) 

2 

16705 

21682 (4.27) 

1 (finest) 

66177 

89550 (4.13) 


The discretization errors on uniformly refined grids at z = 0 and x > 1 are shown in Fig. 62. The finest grid 
denoted as “Grid 1” in the figure corresponds to the level 1 grid in Table 4. The reentrant corner position is at x = 1. 
The same errors on a logarithmic scale are shown versus linearly scaled x (Fig. 62 (a)) and versus the logarithm of 
x — 1 (Fig. 62 (b)); the latter emphasizes errors near the reentrant corner. As expected, the discretization error shows 
first-order convergence at a fixed distance from the corner. The maximum error occurs at the first point away from the 
corner and converges with a \[h order. 

The discretization errors on comer-refinement grids at z = 0 and x > 1 are shown in Fig. 63 and 64 for 2 and 16 
seed cells, respectively. For both sequences of corner-refinement grids, the discretization error shows uniform second- 
order convergence — at all fixed distances from the corner and in the maximum, which again occurs at the first point 
away from the corner. 
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(a) Linear scale in x. 


(b) Logarithmic scale in x — 1. 


Figure 62. Discretization errors in grid refinement at z = 0 and x > 1 for uniform grid. 
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(a) Linear scale in x. 


(b) Logarithmic scale in x — 1. 


Figure 63. Discretization errors in grid refinement at 2; = 0 and x > 1 for corner-refinement grids formed from 2 seed cells. 
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Figure 64. Discretization errors in grid refinement at z — 0 and x > 1 for corner-refinement grids formed from 16 seed cells 


The local discretization errors at x = 1.5 and 2 = 0 for uniform and corner-refinement grids are shown in Fig. 65. 
The logarithm of the errors are shown versus the logarithm of h max in Fig. 65 (a) and versus the logarithm of an 
effective h based on the number of degrees of freedom , h e ff = TV -1 / 2 , in Fig. 65 (b). Triangles indicating the first- 
and second-order convergence slopes are shown for reference. Computations on the uniform grids exhibit first-order 
convergence and computations on corner-refinement grids exhibit second-order convergence on finer grids. The errors 
reduce immediately on both corner-refinement sequences. At fixed degrees of freedom, N, the corner-refinement 
errors with 2 seed cells (corner refinement affects 44% of the domain) are less than the errors with 16 seed cells 
(comer refinement affects 6% of the domain). 




(a) Discretization error convergence versus hmax ■ (b) Discretization error convergence versus h e ff. 


Figure 65. Convergence of discretization errors at a? = 1.5 for uniform and non-uniformly refined grids . 


A parametric study for the corner refinement method has been conducted on a computational domain including 
only the right half of the domain shown in Fig. 61 (a). Dirichlet boundary conditions are prescribed at x = 1,-1 < 
z < 1. The number of recursive cell divisions is varied from two to four; corresponding smallest mesh size obtained 
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after one step of mesh refinement is 1/4 s, l/8s, and l/16s, respectively, where 5 is the size of the coarse comer cell. 
The discretization errors on uniform grids and on grids with corner refinement are shown in Fig. 66. The convergence 
order of the maximum errors (Fig. 66 (a)) improves to the first order with only two recursive cell divisions, but four cell 
divisions are needed to obtain second-order convergence. The L\ norm of discretization errors (Fig. 66 (b)) converges 
with (nearly) second order with two recursive cell divisions. 

The computations indicate that second-order convergence of discretizations errors in all norms can be recovered 
on non-uniformly corner-refined grids with an increased refinement rate in the vicinity of the geometric singularity. 
This refinement is a consistent refinement 36 ensuring that every cell is refined with a rate greater or equal to the rate 
corresponding to the uniform refinement. Although not shown, similar degradation of convergence order to 0[\[Yi) 
locally and to 0(h) globally has been observed on uniformly refined grids for higher-order discretization methods. 32 
To recover design-order convergence for high-order methods on domains with geometric singularities, more local 
subdivisions may be required. To recover design p - th order convergence of the norm of the discretization error, 
the estimated number of subdivisions is 2 p, i.e., each corner cell should be divided into 2p+ 1 cells in each dimension. 
For design-order convergence of the L\ norm, p subdivisions should suffice. Grids used in the corner refinement 
employ more degrees of freedom than corresponding uniformly refined grids, but asymptotically the total number 
of degrees of freedom increases with the same rate as on uniform-refinement grids. The rectilinear grids used for 
the comer refinement place many additional degrees of freedom away from the geometric singularity. These remote 
degrees of freedom seem wasteful; a local adaptive grid refinement is expected to be a much more efficient way to 
improve convergence and to recover the design-order accuracy. 




(a) Loo norm. 


(b) L i norm. 


Figure 66. Discretization error for reentrant corner with uniform grids and grids with varying corner refinements; lj = 2tt. 

The corner-refinement grids used herein are closely related to the optimal meshes developed by Yano and Darmo- 
fal 37,38 for a series of L 2 error control problems. The optimal meshes are derived using established techniques based 
on the polynomial interpolation theory and calculus of variations. For the reentrant-corner problem considered here, 
the optimal meshes are characterized as 

h = Cr k , 

where C is a non-zero constant, f is the distance from the corner, and k is a grading coefficient given by 


p + 1 

The grading coefficient characterizes mesh size disparity. A uniform mesh corresponds to k = 0, and the grading 
becomes stronger (k approaching 1) as a decreases or p increases. Note that k = 0.5 for second-order convergence 
(p = 2) and a zero-angle reentrant corner (a = 0.5). 

AID primal mesh on the interval 0 < x < 1 can be generated with this gradation using the ID mapping 

r = £ 1 ^ , 
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where f = 1 — x is the distance from the corner, and £ is a mapping coordinate, 0 < £ < 1. Discrete meshes are 
obtained using equal spacing in the mapping coordinate, 

Xi = 1 ^ _ 1 1 for i = 1, 2, N, 


where N is the number of nodes on the ID interval. The local average mesh spacing and stretching factor of two 
optimal mesh sequences are shown in Figs. 67 and 68. The gradation parameters k = 0.50 and k = 0.75 correspond 
to p = 2 and p = 5, respectively, for a = 0.5. The average mesh spacing, h avg , and the stretching factor, /?, are 
defined pointwise as 




avg 


)t = 


2 


*^2 + 1 


for i = 2, ..., N — 1. 


The mesh sequences are obtained by doubling the number of cells in each finer grid. Both mesh sequences corre- 
sponding to k = 0.50 and k = 0.75 are consistent, i.e., the mesh spacing at any location goes to zero in the limit as 
N goes to infinity. Also, at a fixed non-zero distance from the corner, the stretching factor goes to unity in the limit as 
N goes to infinity, and the mesh spacing reduces by a factor of two in each refinement. However, the stretching factor 
/?2 corresponding to the interior point nearest to the corner does not approach unity in mesh refinement, being 3 and 
15, respectively, as shown in Figs. 67 and 68. These fixed stretching factors are a consequence of the singularity in the 
mapping at £ = 0. 




(a) Average mesh spacing. 


(b) Stretching ratio. 


Figure 67. Sequence of grids with optimal grading of Yano and Darmofal; k = 0.5; s = 1 — x 

The corner-refinement grids are characterized on average by a grading h = Cr 1_1 / n , where n is the number of 
recursive subdivisions. Similar to the optimal meshes, the corner-refinement mesh sequences are consistent. The grids 
are irregular — the stretching factors vary between one and two in the interior. The corner-refinement grading with 
n = 2 is similar to the grading of the the optimal meshes with k = 0.5 (p = 2); likewise, the corner-refinement 
grading with n = 4 is similar to the grading of the optimal meshes with fc = 0.75(p = 5). 

Particular corner-refinement and optimal grids are compared in Fig. 69 (a) and (b) for grading coefficients of 
k = 0.5 and k = 0.75, respectively, for a comparable number of nodes. The corner-refinement grids start from 2 seed 
cells and have a a portion of the domain that is uniformly refined. The optimal grids are smoother because they are 
mapped grids and provide a slightly lower minimum spacing for a fixed number of nodes N. Although not shown, we 
have verified that the optimal meshes for k = 0.5 recover discretization-error convergence of second order in the L\ 
norm, but only first order in the norm. Computations using optimal meshes with the higher gradation, k = 0.75, 
are expected to recover the second-order convergence in all norms. 
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(a) Average mesh spacing. 


(b) Stretching ratio. 


Figure 68. Sequence of grids with optimal grading of Yano and Darmofal; k = 0.75; s = 1 — x. 


k=0.5 



(a) Corner-refinement grids with n = 2 and optimal grids 
with k = 0.5; 


k=0.75 



(b) Corner-refinement grids with n = 4 and optimal grids 
with k = 0.75. 


Figure 69. Averaged mesh spacing on corner-refinement and optimal grids at equivalent gradation and comparable number of nodes; 

s = 1 — x. 
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